The Hartree-Fock phase diagram of the two-dimensional electron gas 
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We calculate the ground state phase diagram of the homogeneous electron gas in two dimen- 
sions within the Hartree-Fock approximation. At high density, we find stable solutions, where the 
electronic charge and spin density form an incommensurate crystal having more crystal sites than 
electrons, whereas the commensurate Wigner crystal is favored at lower densities, r s > 1.22. Our 
explicit calculations demonstrate that the homogeneous Fermi liquid state - though being an exact 
stationary solution of the Hartree-Fock equations - is never the Hartree-Fock ground state of the 
electron gas. 

PACS numbers: 71.10.-w, 71.10.Ca, 71.10.Hf, 71.30.+h, 03.67.Ac 
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I. INTRODUCTION 

Electrons are found everywhere in matter, most of the time localized by positive charges. In typical condensed 
matter situations, electronic densities and temperatures are such that, in addition to the external positive charges, a 
quantum description of electrons interacting with each other is necessary, leading in general to a difficult quantum 
many-body problem. The homogeneous electron gas, where the positive charges are reduced to solely ensure global 
electro-neutrality, is one of the most fundamental model to study electronic correlation effects. In three dimensions, 
d — 3, valence electrons in alcaline metals realize the electron gas to high precision, in particular in solid sodium 1 , 
whereas the two dimensional electron gas, d = 2, (2DEG), and its extension to quasi- two dimensions^ is relevant 
for electrons at heterostructures, e.g. semiconductor-insulator interfaces^. At zero temperature, the electron gas is 
described by a single parameter, the density n or equivalent ly by the dimensionless parameter r s = a/a#. Here 
a = [2(d — l)irn/ d} -1 ^ is the mean inter particle distance, and a# = H 2 / '(me 2 ) is the Bohr radius, where — e and m 
are the electronic charge and mass, respectively. 

As pointed out by Wigner^, at low densities and zero temperature, electrons will form a crystal, supposed to melt 
at higher densities where the kinetic energy dominates over the interaction. In the limit r s —> 0, the Hartree-Fock 
approximation (HF) applies. Since the non- interacting Fermi sea remains a stationary solution of the Hartree-Fock 
equations, it is natural to assume a Fermi liquid phase at high densities. First principle calculations, such as Quantum 
Monte Carlo^S have located the transition from the Wigner crystal (WC) to the homogeneous Fermi liquid (FL) to 
high precision. Still, there are ind icatio ns that the Fermi liquid phase is not necessarily the absolute ground state of 
the electron gas at high densities ^ 11 * 16 1 and that a direct transition between Wigner crystal and a homogeneous Fermi 
liquid cannot occur in two dimensions in the thermodynamic limit! 12 " 14 i These conjectures actually hold already for 
the electron gas in the Hartree-Fock approximation, but, despite the early predictions by Overhauser of the spin and 
charge density instability of the Fermi liquid ground state, explicit, numerical HF calculations^ have not confirmed 
them for a long time. Based on Bloch functions, these HF calculations^ studied unpolarized and polarized Wigner 
crystal phases of square and triangular symmetries and found a first-order transition to the unpolarized Fermi gas 
which - within this study - remains the lowest energy state for r s < 1.44. Only recently, the first self consistent 
Hartree-Fock solutions with energies below the Fermi liquid energy have been found at high densities 11 . 

The HF solutions of RefP^ obtained without imposing any periodicity in the density show that the fully polarized 
electron gas in two dimensions forms a periodic charge density with triangular symmetry at high densities. In contrast 
to the low-density Wigner crystal, the number of maxima of the charge density is higher than the number of electrons, 
having thus metallic character, and we will refer to such states as incommensurate crystals in the following. However, 
incommensurate states give rise to important size effects, and the calculations in Ref. ^ were limited to ~ 500 
electrons. 

In this paper, we extend the desctiption based on Bloch waves to study arbitrary modulation and occupation 
number. We focus on the density region r s < 4, where incommensurate states may occur. We show how the 
incommensurate states can be represented by the vector Q of the charge modulation. Restricting the search for 
the HF ground state to states with arbitrary Q, we are able to overcome size restrictions and we explore the phase 
diagram of the 2DEG including triangular (A) and square (■) symmetries. While our minimization also includes 
the possibility of partial polarized states, they do not occur as ground states which are either unpolarized (U) or 
fully polarized (P); in particular, we show that the incommensurate unpolarized crystal is favored at high densities. 
Whereas the momentum distribution of the Wigner crystal is a continuous function of the momentum, we show that 
there are angle-selective steps in the incommensurate phase. 



II. METHODS 

The Hamiltonian of the electron gas containing N electrons writes 

H = --J2^i+ Yl v ( x i~ x j) (!) 

i l<i<j<N p 

where is the Laplacian with respect to v(x) is the electrostatic interaction v(x) = and we have used 

atomic units where distances are measured in units of a# and energies in Hartree, lHa = h 2 / \ma 2 B ) . In addition 
to Eq. 0, the interaction between electrons and a positive background charge must be considered to ensure charge 
neutrality. 

We are considering N electrons in a finite box of volume V, of sizes L\ and L2, with periodic boundary conditions, 
so that the momentum k belongs to the lattice L* generated by L\ and L\ satisfying LiL* = 27r5^ . 
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FIG. 1. (color on line) Illustration of the /c-space in the square (left) and triangular (right) geometry. At the center of each 
figure, are shown the brillouin zone B (in white) and the corresponding basis vectors Qi and Q2. The first and second shell 
of neighboring cells, B + n\Q\ + 712Q2, are shaded in light-green and light blue, respectively. For square (resp. triangular) 
symmetry, the integers m of the first and second shell satisfy n\ + n% — 1,2 (resp. n\ + 7t| — 71^2 = 1,3). The corresponding 
number of cells are summarized by M\ in the right column. Most of the results presented in this paper are done including 
a number of bands, M A , which corresponds to two neighboring shells for the square and one for the triangular geometry. In 
light-red, we indicate the elementary cell, Bo, used in our numerical calculations; the N B black dots are an example of the 
discretization of the Brillouin zone (here N B = M 2 with M — 8, as explained in Sec III). The circle indicates the Fermi surface 
of a Fermi gas. 



Within the Hartree-Fock approximation, the energy expectation value is minimized with respect to a single skew- 
symmetric product of N single particle states. Periodic solutions are special states which can be described by Bloch 
waves. Let A be a sub-lattice of L* generated by Qi and Q 2 . The Brillouin zone is defined as the Voronoi cell of the 
origin and a periodic state is given by \(fk) = X^ga a k(q)\k + <?), where k belongs to the Brillouin zone, B (see Fig{l]). 

As a particular case, the Wigner crystal (WC) is obtained by choosing A such that the Brillouin zone B contains 
exactly N states where N is the number of electrons. Thereafter, the state is built as f\keBWk)- An upper bound 
of the ground state energy is obtained by minimizing the coefficients a^iq) of the Bloch functions. As r s approaches 
zero, the kinetic energy dominates which is minimized by the Fermi gas (FG) defined by k < kp (see below). Such 
a state cannot be described by the WC. In refP^, the HF energy of the polarized gas has been minimized without 
imposing periodicity of the solutions; nevertheless, at intermediate densities, the HF ground states are periodic with 
larger modulations than in the WC corresponding to less electrons than states in the Brillouin zone. 

In this paper we focus on periodic solutions with arbitrary modulations. For a given modulation and for fixed choice 
of /c-vectors in the Brillouin zone the energy computation is fast enough to tackle millions of electrons as every single 
state is described by very few parameters only. However, the minimization with respect to the choice of /c-vectors in 
the Brillouin zone becomes a complicated combinatorial problem. This combinatory problem is simplified within the 
framework of density matrix. 

The one body density matrix, pi, is a symmetric positive matrix such that Trpi = 1. Provided that pi < 1/iV, this 
matrix can be seen as one body density matrix of a state of TV electrons. In the thermodynamic limit, the two body 
uncorrelated density matrix can be expressed in terms of pi as 

p 2 (l, 2; 1', 2') = pi(l- 10/91(2; 2') - Pl (l; 2>i(2; 1'). (2) 

The total energy, a-priori a function of the reduced one and two body density matrices, can be expressed entirely as 
a functional of pi. Explicitly, we obtain for the energy per particle in atomic units 

E =1, ^2 k 2 pi(k,(T',k,(T) + ^ ^2 v q p 2 (k 1 ,a 1 ,k2,02\ki+q,G 1 ,k 2 -q,(T2) (3) 

ci ,cr 2 

where v q = l/\\q\\ for q ^ 0, and = 0- For instance, the unpolarized Fermi gas (U-FG) corresponds to pi(ka, k'cj') = 
Skk'S**>Q(kF,u ~ \\k\\)/N with nk 2 FU = 2ir 2 N/V; the resulting energy is E U FG = l/(2r 2 s ) - 8/(37r>/2r s ), and the fully 
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polarized Fermi gas (P-FG) corresponds to pi(ka, k f a f ) = Skk , S a +S a(T fQ(kF,p — \\k\\)/N with i\k\ p = (2tt) 2 N/V and 
energy E PG = 1/r 2 . — 8/(37rr s ). In general, without any specification, kp denotes the Fermi wave vector according to 
the polarization of the corresponding state. 

In the following we restrict the density matrix to represent periodic solutions. The corresponding one-body density 
matrix can be written as : 

Pi(M') = Pi(k + q,a;k + q',a f ) = p k (q, cr; q' , cr') (4) 

with g, q' G A and k in the Brillouin zone B. Thus, the density matrix is now described by a family of positive matrices 
Pk such that pk < l/N and ^ k Tr pk = 1. 

Numerically, we truncate the number of lattice vectors of the sub-lattice A and include only the first M A vectors of 
smallest norm in the numerical calculations. In the framework of band structure calculations, where the Bloch states 
are obtained from an external periodic potential, M A corresponds to the number of bands considered. Thus pk is a 
2M A x 2M A matrix and in order to fulfill the condition pk < 1/iV, it is more convenient to write: 

p k = U^D k Uk (5) 

where Dk is a diagonal matrix with < Dk < l/N and Uk is a unitary matrix. The potential energy contains a 
convolution in momentum space calculated using fast Fourier transform (FFT). The minimization of the HF energy 
is done computing the gradient of the energy with respect to Uk and Dk. The only drawback of the method is to 
fulfill the condition D k < l/N. 

The minimization at given density consists in the following steps. At first we choose Dk and Uk to start with. Then 
we find the best Uk with a quadratic descent methocP4 The next step is to try to improve Dk given the gradient 
of the energy with respect to Dk and the linear constrains, < Dk < l/N and ^2 k Dk = 1. The process stops as 

soon as D ( ™ w) =D k .ln this case almost every Dk are or l/N and the gradient is negative or positive accordingly. 

Otherwise, we change Dk into (1 — e)Dk + eD^ ew ^ (with a small e to ensure that Uk follows Dk adiabatically) and 
we restart the minimization with respect to Uk- 

In this work, we study the 2DEG for triangular (a) and square (■) symmetries where ||Qi|| = HQ2II = Q- Starting 
from a state of arbitrary polarization, the minimization always resulted in either an unpolarized (U) or a fully polarized 

(P) state. The Brillouin zone of the Wigner crystal contains exactly N states, so that Q/kp = Qw /kp = \J 2tt/^/3 « 
1.9046 for the triangular WC (U or P), whereas Q/kp = Q w /k F = y/n « 1.7725 for the square WC (U or P). 
Notice that for triangular symmetry the corresponding direct space lattices are quite different: honeycomb lattice for 
unpolarized and triangular lattice for polarized states. The FG can be reached when the Fermi surface is contained 
inside the Brillouin zone, that is for Q > 2kp. Thus, in our simulations, Q varies between Qw and 2kp. 

III. CONVERGENCE STUDIES 

We first focus on size effects in the thermodynamic limit extrapolation, N — >• 00. We set Qi — ML*, thus the 
Brillouin zone contains N B = M 2 vectors. Since N/N B = (Qw/Q) 2 , this limit at fixed Q is equivalent to study the 
convergence with respect to N B . Fig. [2] shows the size extrapolation of the 2DEG in the triangular symmetry at 
r s = 4 (Q = Qw), together with the results of Trail et alP^, done at N B = 13 and M A ~ 20, and those of Ref.^4 As 
the calculations of Ref. 11 do not assume any periodicity in the HF search, they are limited to system sizes TV < 500, 
and the extrapolation to the thermodynamic limit is less accurate. 

Size effects depend on the phase considered. In the incommensurate phase, size corrections are not any more 
monotonic functions, as in the Wigner crystal, but oscillatory behavior occurs depending on the density r s , and on 
the modulation vector Q. In Figj3j we show the energy of the 2DEG in a triangular symmetry at r s = 2.5 versus 
the modulation Q (incommensurate crystal) for various system sizes using M = 2 P , with p from 4 to 9 (N B = 16 2 
up to 512 2 ). Note the random like oscillations due to the discretization, N B , of the Brillouin zone. However, at large 
enough N B , these oscillations are sufficiently small to analyze safely E(Q,r s ) as seen in FigjE] 

Our second parameter is the number of vectors M A considered in A. Note that truncation of A does not violate 
the variational principle, so that the energy of a converged HF solution must decrease as M A increases. Figure [2] and 
Fig. [4] show the convergence in system size N B (discretization of the Brillouin zone) together with the exponential 
convergence in M A which measures the large k importance. As expected, energies decrease with M A because the 
Hilbert space is increased. Interestingly, the M A improvement is mainly independent of N B (see Figj2]-right and 
Figj4|, which allows us to work with small M A and estimate corrections using small systems. Most of the calculations 
presented in this paper are thus performed with M A = 7 and M A = 9 bands for supercells of triangular and square 
symmetry, respectively. 
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FIG. 2. (color on line) Energy (in Hartree units) of the 2DEG in the triangular symmetry at r s — 4 (WC) as a function of 
the number of particles, N = N B , and the number of included bands, M A . Em = — 1-1061 / r s is th e Madelung energy. (PA 
indicates polarized final state with triangular geometry.) Left: comparison with previous worl J 11 * 15 ! 17 ! Blue full down triangles 
are results of the present work using M A = 19. Right: convergence with respect to N and M A . The inset is a zoom of the 
dotted- line-domain. 




IV. RESULTS 



We have studied the HF ground state of the 2DEG in the density region 0.8 < r s < 4 at zero temperature 
considering commensurate and incommensurate solutions with square and triangular symmetries. At low densities 
the electro ns form a commensurate Wigner crystal of modulations Q = Qw and we recover the results of previous 
HF studies ^ * 15 * 17 *. For higher densities, an incommensurate crystal with modulation Qw < Q < 2kp is formed for 
any fixed polarization and symmetry. 

Figure [5] summarizes the energy gain with respect to the unmodulated Fermi gas as a function of Q at different 
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FIG. 4. (color on line) Convergence of the energy with respect to M A for the 2DEG with triangular symmetry for two system 
sizes M = 32 and M = 64. Full and dotted lines stands AE = Em a x — Em a 2 an d AE = Em a 2 — Em a 3 , respectively, with 
Ma,i = 7, Ma,2 = 13, Ma,3 = 19. Crosses and dots stand for M = 32 and 64, respectively. (Values at r s = 2 are close to the 
convergence threshold of the descent method.) 




FIG. 5. (color on line) Energy difference, with respect to the Fermi gas, E(Q,r s ) — ^g P , in milli Hartree versus modulation, 
Q, for different densities and symmetries (A or ■) at N B — 256 x 256. The final polarization obtained after minimization is 
either unpolarized (U) or fully polarized (P). Lines are the polynomial fits using the parameters given in Table-]!] In each figure, 
the lowest curve (largest r s ) with triangular or square symbols has a minimum at Q = Qw- The bold dashed-line connects 
Qm(v s ), the minima of E(Q,r s ) for fixed r s . Vertical dotted lines indicate Qw- 



densities. Well inside the incommensurate phase (Q > Qw), the energies can be well represented with a polynomial 
form: 

3 2 

E(Q, r s ) = E FG {r s ) + ^ ^ ^ XV * ( 6 ) 

i=0 j=0 

where X = 100(Q/kp — 2). The parameters ol^ determined by least square fits are given in Table-|l| From this 
parametrization, for fixed r s , we determine the minimum Qm^s) of E(Q,r s ), shown in FigjHJ 

The incommensurate phase is characterized by a crystal in direct space with slightly more lattice sites N B than 
electrons N, increasing for larger modulation according to N B /N = (Q/Qw) 2 - Figure |6] shows typical charge and spin 
densities in the incommensurate phase for the triangular and square geometry. The two examples are chosen close 
to the transition to the Wigner crystallization. The amplitude of the modulation of the charge densities is about an 
order of magnitude smaller than that of the spin densities, an effect which is even more pronounced at higher density. 

The momentum distribution rtk (N times the diagonal part of pi) provides additional insight. In contrast to 




FIG. 6. (color on line) One body charge and spin densities of an unpolarized incommensurate crystal with triangular symmetry 
(left: r s = 1.2, Q/k F = 1.933, N/N B ~ 0.97) and square symmetry (right: r s = 1.5, Q/k F = 1.844, N/N B ~ 0.92). Average 
values have been subtracted. Lengths are given in units of the inverse modulation, Q -1 . The color scaling is the same for all 
pictures. Contour levels are ±0.01, ±0.02 for the charge densities and at ±0.1, ±0.2 for the spin densities . 
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FIG. 7. (color on line) Modulus of the momentum distribution |n^(fc)| = |n^(/c)| of U solutions as a function of wave vector k 
in the positive quadrant (other parts can be deduced by symmetry). Top line and bottom line are for triangular and square 
symmetry, respectively. Contour levels are at 0.5 and 0.1, 0.01, etc. From right to left is shown the evolution from the Wigner 
crystal distribution (continuous function everywhere) to Fermi gas with a step along some directions. For the incommensurate 
states, note the step-function behavior to a domain where nk — which grows in the corner of the Brillouin zone when r s 
decreases. 



the step-function behavior at kp of the Fermi gas, is continuous inside the commensurate Wigner crystal phase 
and its variation reflects the symmetry of the Brillouin zone. The incommensurate phase still reflects the underlying 
symmetry of the crystal, but angle selective steps occur at the corners of the Brillouin zone (see Fig. [7]). The rounding 
of the corners increases for smaller r s , and the isotropic step-function of the Fermi gas is continuously approached for 
rwO. 

Whereas we have found that the incommensurate phase is always favored compared to the Fermi gas solution, 
independently of the imposed polarization and crystal symmetry, the unpolarized incommensurate hexagonal crystal 
becomes the true HF ground state at high densities, r s <r c s ~ 1.22. The different phases and energies for 0.8 < r s < 4.0 
are illustrated in FigjH] Although our HF method does not impose the polarization, we have not found any stable 
partially polarized ground states. At r s > r c s the unpolarized electrons form a commensurate Wigner crystal of 
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FIG. 8. (color online), (a): Phase Diagram of the 2DEG at T = 0, where Em = — 1.1061/r s is the Madelung energy and 
energies are multiplied by r Z J 2 • Dotted lines correspond to the Fermi gas. Blue and red curves represent the triangular (A) 
and square (■) phases, and P and U stands for polarized and unpolarized phases, respectively, whereas W and I (thick curves) 
indicate Wigner crystal and incommensurate crystal, respectively. The dashed vertical lines indicate the transitions. (b): Energy 
gain with respect to the unpolarized Fermi gas energy, Ep G , in the high density region. 



TABLE I. Coefficients ciij of the polynomial fits E(Q,r s ) — Efg(t s ) defined by Eq[6] 



U-A 


U-B 


P-A 


-0.78611 1.88240 -1.13180 
0.35614 -0.64435 0.34780 
0.10624 -0.15804 0.05531 

-0.00166 -0.00044 0.00081 


-0.621900 1.53040 -0.96941 
0.058858 -0.26652 0.27359 
0.032321 -0.06425 0.02500 

-0.022986 0.02665 -0.00750 


0.15758 -0.08577 -0.011495 
0.24875 -0.24028 0.070520 
0.11155 -0.10668 0.024822 
-0.00090 -0.00321 0.001310 



hexagonal symmetry, and, at r s ~ 1.62 a structural transition from the unpolarized hexagonal WC to the unpolarized 
square WC occurs, followed by a transition from the unpolarized square WC to the fully polarized triangular WC at 
r s ~ 2.6. 



V. CONCLUSION 



We have studied the 2DEG in the Hartree-Fock approximation at densities r s < 4. We confirm previous observations 
of incommensurate phases of the fully polarized electron gas^, performing calculations of much larger system sizes. 
We further included electron polarization, as well as square and triangular symmetries. Our HF phase diagram at zero 
temperature is much richer than that obtained previously^, which did not consider the unpolarized triangular WC, 
nor any incommensurate phase. Our numerical calculations explicitly confirm the old conjecture of Overhausei^^ 
that Fermi gas is never the HF ground state which has been proven rigorously for the fully polarized electron gadm 

We have further shown that the momentum distribution provides an unambiguous characterization of the incom- 
mensu rate p hase. In contrast to the isotropic momentum distribution of a Fermi liquid, discontinuous at the Fermi 
surfac e 18 * 19 !, the incommensurate phase exhibits an anisotropic momentum distribution intermediate between a crystal 
and the Fermi gas with forbidden domains inside the Brillouin zone, where rik jumps to zero. 
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